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ABSTRACT 

The study of warm molecular gas in the inner regions of protoplanetary disks is of key importance 
for the study of planet formation and especially for the transport of H2O and organic molecules to 
the surfaces of rocky planets/satellites. Recent Spitzer observations have shown that the mid- infrared 
spectra of protoplanetary disks are covered in emission lines due to water and other molecules. Here, 
we present a non-LTE 2D radiative transfer model of water lines in the 10-36 fim range that can be 
used to constrain the abundance structure of water vapor, given an observed spectrum, and show that 
an assumption of local thermodynamic equilibrium (LTE) does not accurately estimate the physical 
conditions of the water vapor emission zones, including temperatures and abundance structures. By 
applying the model to published Spitzer spectra we find that: 1) most water lines are subthcrmally 
excited, 2) the gas-to-dust ratio must be as much as one to two orders of magnitude higher than the 
canonical interstellar medium ratio of 100-200, and 3) the gas temperature must be significantly higher 
than the dust temperature, in agreement with detailed heating/cooling models, and 4) the water vapor 
abundance in the disk surface must be significantly truncated beyond ~ 1 AU. A low efficiency of 
water formation below ~ 300 K may naturally result in a lower water abundance beyond a certain 
radius. However, we find that chemistry, although not necessarily ruled out, may not be sufficient to 
produce a sharp abundance drop of many orders of magnitude and speculate that the depletion may 
also be caused by vertical turbulent diffusion of water vapor from the superheated surface to regions 
below the snow line, where the water can freeze out and be transported to the midplane as part of 
the general dust settling. Such a vertical cold finger effect is likely to be efficient due to the lack of a 
replenishment mechanism of large, water-ice coated dust grains to the disk surface. 
Subject headings: astrochemistry - line: formation - planetary systems: protoplanetary disks - radia- 
tive transfer 



1. INTRODUCTION 

Planets are thought to predominantly form within 
10 AU of young low mass stars (< a few solar masses). 
Concurrent with the process of planet formation, a rich 
carbon, nitrogen and oxygen chemistry evolves to even- 
tually form the building blocks of the icy bodies within 
young planetary systems; moons around giant planets, 
comets and Kuiper belt objects. This chemistry will 
also seed the surfaces of rocky terrestrial planets with 
water and organics, probably either through impacts 
with comets or carried in hydrated refractory meteorites 
(|Ravmond et al.ll2004 ). 

This early, highly active and evolving stage takes place 
within the first few million years of a star's lifetime while 
the protoplanetary disk is still gas-rich. The planet form- 
ing region can therefore be traced by observations of 
warm (100-2000 K) molecular gas, through a multitude 
of rotational and rovibrational transitions spanning the 
infrared wavelength range (2-200 /zm). 

One critical question is how the distribution of wa- 
ter occurred in the Solar Nebula and how th e surfa ce of 
the young Earth was hydrated. ISalvk et all (|2008l ) and 



ICarr fc Najital (|2008f l report detections of hundreds of 
water vapor lines in three protoplanetary disks observed 
with the Spitzer Space Telescope InfraRed Spectrome- 
ter (IRS) between 10 and 36 /jm. These first, pioneering 
spectra were analyzed using simple 'slab models' that as- 
sume level populations consistent with thermodynamic 
equilibrium and otherwise consisting only of a column 
density, a single excitation temperature and an emitting 
solid angle. 

Fits to data using such simple slab models result in 
column densities of Nco = 6 — 7 x 10 18 cm~ 2 , Ah 2 o = 
6 - 8 x 10 17 cm" 2 , and A OH = 2 x 10 17 cm -2 and an 
excitation temperature T — 1000 K for DR Tau and AS 
205A, which is simila r to the observations of AA Tau by 
ICarr fc Najital (|2008h . The latter paper also reports the 
detection of CO2 and that of the small organics C2H2 and 
HCN. Given a slab model, the line emission is derived to 
be confined within radii ranging from 0.6 to 3.0 AU. 

Recently, a number of theoretical models have 
appeared that focus on the gas-phase thermo- 
chemical modeli ng of the inn er re gions of proto- 



plane t ary disks dMarkwick et al.l l2002t iG lassgol cfet alj 
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adopted physics (e.g. 
FUV, the possibility of transport throughout the disk, 
time-dependent treatment of chemistry) and deal with 
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various chemical aspects (e.g., organic species, H2O, 
carbon fractionation) . Other theoretical models focus on 
the dynami cal transport of ind i vidual species through 
the disk (IStevenson fc Lunind Il988t ICiesla fc Cuzzil 
200$ ICieslal 12009ft . A kev question is therefore whether 
the large amounts of warm (T ~ 300 — 1000 K) water 
that have been observed can be explained by models 
with either in situ formation or dynamical transport 
from larger radii. 

To date, neither the observational work nor the the- 
oretical papers that consider the chemistry in the inner 
region of protoplanetary disks have shown how physical 
parameters - specifically the actual abundance structure 
of a given molecule - can be derived from line observa- 
tions in a meaningful way. In this paper, we propose 
that, given the complexity of inner disk chemistry and 
volatile dynamics, progress can best be made by sepa- 
rating the chemical-dynamical model from the radiative 
transfer modeling of observations. This way, high qual- 
ity observations, such as infrared imaging spectroscopy, 
can be used to derive an abundance structure of a given 
molecule by application of a radiative transfer code that 
performs: 1) a non-LTE excitation calculation in a disk- 
geometry followed by, 2) raytracing that correctly sam- 
ples the velocity field at small radii, and 3) a front 
end that simulates the action of the appropriate tele- 
scope/instrument. Given a derived abundance structure, 
the interpretation step, using chemical-dynamical mod- 
els, becomes much easier than the approach in which one 
begins with the chemical model. This paper describes 
step 1). A companion p aper describes steps 2) and 3) 
(jPontoppidan et al.ll200g) . 

Here our major focus is an examination of the excita- 
tion and infrared radiative transfer of water vapor in the 
surface layers of protoplanetary disk. We find that water 
vapor must be strongly depleted from the surface over a 
significant range of radii. This can be due to a chemi- 
cal effect, as wate r is less efficiently form ed at tempera- 
tures T < 300 K (jGlasseold et al.l l2009D ; but we specu- 
late that additional depletion is active due to vertical gas 
diffusion, followed by freeze-out below the snow line and 
settling to the mid-plane on icy dust grains. This pro- 
ce ss is analogous to the r adial cold finger effect proposed 
bv IStevenson fc Lunind (|1988l ). but will not be balanced 
by th e inward radial migration of icy bodies (|Salvk et al.l 
2008). To distinguish it from the radial cold finger effect, 
we call it the vertical cold finger effect (see Section 15. 2|) . 

To place the non-LTE model into context, we first de- 
scribe slab models of disk infrared emission and the col- 
lisional/radiative excitation of water. The details of the 
statistical equilibrium calculations are then presented, 
along with our initial results and general conclusions. 
The models presented here are meant only to provide a 
fundamental analysis of the global properties of water 
vapor emission from disks, the more elaborate parame- 
ter study of the water excitation model that is needed to 
fit the Spitzer IRS spectra of individual sources is post- 
poned to a separate paper (Meijerink et al. 2009, in 
preparation) . 

2. LTE SLAB MODELS VERSUS NON-LTE 2D MODELS 

Recent pap ers from ISalvk et "ail (|2008l ) and 

ICarr fe Naiital (|2008l ) have been successful in ob- 
taining good matches to H2O line observations in the 



Mid-Infrared (MIR) Spitzer-IRS band (A ~ 10 — 36^m) 
using single temperature, single column density LTE 
models with no global kinematics, such as Keplerian 
motion. These fits have been used to provide numerical 
estimates of the warm water column densities. Since 
a key point of this work is to show why non-LTE 
calculations are required to quantitatively interpret 
the observed molecular line spectra of the warm inner 
regions of protoplanetary disks, particularly the emission 
from H2O, we first examine the successes and limitations 
of LTE 'slab' models. We will argue below that taking 
the properties derived from such models at face value is 
dangerous and may well lead to erroneous conclusions 
when attempting to compare abundances with chemical 
models. Specifically, a good fit in a given wavelength 
range does not imply that the fitting model is correct. 
The low resolution of the Spitzer spectra, and the result- 
ing lack of line profile information, introduce significant 
degeneracies, as we discuss next. We will then show 
that the introduction of relatively simple physics lead 
to constraints not available from slab models, and that 
a requirement to match all lines in the observed Spitzer 
range can significantly constrain the results. Although 
the lineshapes also contain a great deal of info, it is 
not considered here, since the IRS has a resolution of 
R = A/AA = 600 (500 km/s). 

Slab models ignore the complex geometric structure of 
protoplanetary disks. In the inner regions, densities can 
be as high as tih > 10 16 cm~ 3 in the midplane and at 
small radii (R < 0.5 AU), but decrease with radius and 
height (see Fig. [J) over many orders of magnitude to den- 
sities n < 10 3 cm" 3 . Gas temperatures can be as high 
as T ~ 5000 K in the u pper reaches of the d isk atmo- 
sphere where the X-rav (IGlassgold et al.ll2004h or FUV 
(|Kamp fc Dullemondl l2004) attenuation is minimal. The 
temperature drops steadily in the shielded regions of the 
disk to values of T ~ 10 — 100 K, where the gas tempera- 
ture couples strongly to the dust temperature. Therefore, 
the disk consists of an ensemble of densities and temper- 
atures. Lines are, in general, the result of emission from 
a range of disk radii, and thus cannot be represented by 
a single excitation temperature. The mid-infrared wave- 
length region contains radiative transitions from upper 
levels with an enormous range of energies, from values 
at/near the photon energy (E < 720 K or 500 cm -1 ) to 
values in excess of E > 7200 K or 5000 cm -1 . The con- 
version between wavenumber and temperature units is 
1.44. Temperature units are used throughout this paper. 
For a light, asymmetric top molecule, such as H2O, tran- 
sitions with various upper state energies can be highly 
scattered throughout the spectrum, and even a small 
spectral region can contain lines that trace very different 
regions of the disk. Therefore, the assumption that all 
lines are formed in the same region of the disk is imme- 
diately seen to be inappropriate - and can easily lead to 
an erroneous determination of the line optical depth in a 
slab model. In fact, it is essentially meaningless to even 
assign an optical depth to a given line, since emission 
from different disk radii will have local optical depths 
varying by orders of magnitude. 

Furthermore, LTE assumes that the level populations 
are in accordance with the Maxwell-Boltzmann distribu- 
tion at a particular temperature. This only holds when 
collisions dominate the molecular excitation/relaxation, 
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or numerically when the ambient densities are higher 
than the critical densities n cr it = A u i/C u i, where A u i 
and C u i are the radiative and collisional decay rates of a 
particular transition. The critical densities for the exci- 
tation of the water lines through collisions range between 
ncrit = 10 8 to > 10 12 cm -3 , and tend to increase strongly 
with excitation energy. Pure rotational transitions with 
low excitation energies, which lie in the frequency range 
of the Heterodyne Instrument for the Far-Infrared (HIFI) 
on the Herschel Space Observatory (Herschel), have crit- 
ical densities in the range 10 8 — 10 9 cm~ 2 . Unfortunately, 
these lines will likely be difficult t o observe toward disk s 
with HIFI for sensitivity reasons (|Meijerink et al.l l2008). 
Spitzer-IRS (10 — 36 /mi range) data, on the contrary, 
have shown rich spectra of warm water in the inner re- 
gions of protoplanetary disks. These consist of rotational 
transitions with even larger critical densities than the 
HIFI lines and are therefore unlikely to be thermalized. 
Lower excitation lines of a given species tend to have 
lower critical densities, and as a result, lines with lower 
excitation energies will trace larger regions of the disk. It 
is therefore expected that the line shapes and line images 
will differ significantly from one transition to another. 
Thus, future spectrally and spatially resolved data will be 
crucial for constraining the water abundance structure, 
and should be interpreted with tools that specifically 
model (or mimic) the physical and chemical properties 
of the disk. 

3. EXCITATION OF H 2 

There are a number of ways to excite water to higher 
energy states: 1) collisions with, e.g., H, H 2 , He, and elec- 
trons; 2) excitation by radiation produced by hot dust 
located at the inner rim of the disk; 3) photo-desorption 
from dust grains; and 4) chemical formation in highly ex- 
cited states. In this paper, we only take into account 1) 
and 2). For this study all hydrogen is assumed to be in 
molecular form, and He excitation is considered by scal- 
ing the molecular hydrogen rates. Electron excitation is 
neglected, since electrons arc expected to be abundant 
only in unshielded low density regions, where water is 
subthermally excited and its fractional abundance very 
low, which therefore do not contribute much to the over- 
all line emission. 

Water lines in the 10 — 36 /im wavelength range cov- 
ered by the Spitzer-IRS Short-High (SH) and Long- 
High (LH) modules are overwhelmingly dominated by 
pure rotational transitions between levels in the ground 
(u 1 w 2 W3)=(000) vibrational state as well as between those 
in the first vibrational state (uii;2U3)=(010). Further, for 
a non-LTE calculation, transitions that connect the two 
vibrational levels must also be considered. 

Given the widespread importance of water, there have 
been significant efforts to calculate the H2O excitation 
rates - both experimentally and theoretically. Unfor- 
tunately, most studies only give state-to-state rate co- 
efficients below the fir st excited vibrational level, i.e, 
below E = 2294.4 KjGreenet "all fl993: Phillip s eTall 
19961 : iFaure et alJ I2004L l2007t iD ubernet & Grosieanl 
200llGrosiean et al.ll2003h . The only H 2 0-H 2 and H 2 0- 
electron collisional rates for higher levels of ex c itatio n 
currently available are from IFaure &: JosselirJ (2008). 
This database contains all collisional rotational and ro- 
vibrational transitions between the 5 lower vibrational 



levels, («iu 2 « 3 )=(000), (010), (020), (100), and (001), up 
to an energy E = 7200 K. The authors expect the data 
to be accurate within a factor of ~ 5 for the largest rates 
(> 10~ n cm 3 s _1 ). Smaller rates, which also include 
the rovibrational excitation/de-excitation, are expected 
to be accurate to an order of magnitude. Summed rates 
and cooling rates have better precision and the data are 
well suited for the modeling of emission spectra. 

In Fig. [TJ the intensities for lines at T = 700 K are 
shown, selected to be in the Spitzer-IRS spectral range 
between 10 and 36 /im. The top panel s hows the inten- 
sity o f all lines in the HITRAN database (Rothma n et al.l 
2005), where the lines with upper level energy E < 
7194 K are the open green circles (n = 1454), and the 
lines with E > 7194 K in red filled circles (n = 142). 
Because collisional rates between levels with energy E > 
7194 K are not available it is possible that potentially 
strong lines are not calculated, especially in the 10 - 
15 /im spectral range. Therefore, the calculation of these 
rates is important for future modeling. 

To limit the computational load, the non-LTE calcu- 
lation is restricted to transitions in the (uiu 2 i>3)=(000) 
and (010) vibrational levels. In the lower panel in Fig. 
[1] we show the intensity of transitions with upper level 
energy E < 7200 K, where the transitions between levels 
in the (it ^2^3) =(000) and (010) are in green open circles 
(n = 1230) and all others (n = 224) are in red filled cir- 
cles. It can be seen that the excluded vibrational levels 
contain lines that are weaker by 1-2 orders of magnitude 
relative to lines in the vibrational levels modeled. At the 
Spitzer-IRS SH+LH spectral resolution of R=600, these 
lines are unlikely to contribute significantly to the overall 
spectrum. 

4. MODEL SETUP 

Next we investigate how the H 2 10 - 36 /im spec- 
trum depends on the properties of the disk. The aim 
is to derive a fiducial model that qualitatively matches 
the observed data. Quantitative fits to the spectra of 
individual sources are postponed to a later paper. The 
central star is chosen to be a representative T Tauri. It 
has mass M = 1.0 M©, age t = 2 Myrs, and solar metal- 
licity Z = 0.02, corresponding to a radius R — 2.0 R Q 
and effective temperature Teg = 4275 K whe n using the 
pre-mam sequence tracks of ISiess et all (|2 QQ0j). The ste l- 
lar spectrum is a Kurucz stellar model (jKurucd ll993). 
An inclination of 30° is adopted for rendering of spectra. 

The disk has a mass of M = 0.01 M , outer radius 
i?out — 120 AU, flaring parameter H/R cx i? 2 / 21 , and 
outer pressure scale height h p /R = 0.10. These param- 
eters determine how much light is intercepted at each 
radius and thus the temperature distribution of the disk. 
We adopt a Gaussian density distribution in the vertical 
direction and a radial surface density S cx R . The disk 
is roughly in hydrostatic eq uilibrium (for an isotherma l 
disk), but flatter than the [Chiang fc Goldreichl (|1997l ) 
disk (H/R cx R 2 ' 7 ), giving us a better correspondence to 
the observed spectral energy distributions in our sample 
of T Tauri disks. However, we emphasize that the model 
is not an attempt to fit any particular disk. The dust 
sublimation temperature is set at T ~ 1600 K, giving a 
radius of the inner rim of R r i m = 0.075 AU. 

In addition to parameters determining the dust struc- 
ture of the disk, a number of parameters are needed 
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TABLE 1 
Model input 
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Wavelength [fim] 
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Fig. 1.— LTE intensities at T = 700 K for: (top) Lines with 
upper level energy E < 7194 K (open green circles) and upper 
level energy E > 7194 K (filled red circles); (bottom) Lines with 
lower and upper level in v=000 or v=010 (open blue circles) and all 
other transition (orange filled circles) both with upper level energy 
E < 7194 K. 



for the line modeling. These are the prescription of 
the intrinsic line width, the freeze-out model, the frac- 
tional H 2 abundance and the gas temperature and den- 
sity structure. We assume that the intrinsic lincwidth 
throughout the grid is determined by thermal broaden- 
ing (vth = \/2fcT//iH20 TO p) 5 an d turbulent broadening, 
assumed to be a fraction of the sound speed Vturb = ec s , 
where e = 0.03 is adopted in this paper. The intrinsic 
linewidth is therefore dominated by thermal broadening. 
Five different model types are considered: 1) Constant 
water abundance throughout the disk; 2) Reduced water 
abundance in regions where water freezes out (see Fig. 
[3]); 3) Decoupling of gas temperature and destruction of 
water where gas is exposed to direct stellar irradiation; 
4) Increasing the gas-to-dust ratio from 1280 (which is al- 
ready and order of magnitude beyond the standard inter- 
stellar medium value) to 12800, and 5) Truncation of the 
water abundance beyond R ~ 1 AU due to the vertical 
cold finger effect. The model parameters are summarized 
in Tabled! 

The water emission spectrum is calculated as follows: 
(i) First, the dust temperature distribution and 
the mean intensity in the disk is determined us- 
ing the 2D dust radiation transfer code RADMC 
(Dull emond fc Dominikl I2004D . The mean intensity is 
used to calculate the radiative (de-)excitation rates £? y < 
Jjj >. The adopted dust mixture contains 15 percent 
carbon grains and 85 percent silicate grains. The two 
dust components are thermally coupled. The dust mass 
distribution and opacities are shown in Fig. [51 The 
silicate dust size distribution has a peak at 1 — 10 /im 
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a = 2/21 

0.10 

1.0 

2.0 

4275 

see Fig. |2] 

1280 (1-3); 12800 (4-5) 
30° 



grains, while the carbon grains consists of both small and 
very large grains (> 10 fim). This corresponds to grains 
that hav e grown beyond those foun d in the interstellar 
medium (jKessler-Silacci et al.ll2006| ). It is implicitly as- 
sumed that grains larger than 40 fj,m have settled to the 
midplanc where they do not affect the water lines (which 
are formed near the disk surface). 

(ii) The level populations a re calculated with the ra- 
diativ e excitation code /33D (|Poelman &; Spaani 120051 
2006), a mult i -zone escape probability excitation code. It 
can be applied to arbitrary geometries and is suitable for 
any atom and molecule. This code is adopted because it 
is about 10-100 times faster than existing MC/ALI codes 
- especially when the lines are highly optically thick and 
IR pumped. While /33D can be used in 3 dimensional 
structures, it is only feasible in applications where a lim- 
ited amount of levels (~10) are used. Here we adopt a 
1+1D approximation to alleviate the computational ex- 
pense of iterating over a total of ~ 540 levels for both 
ortho and para water, especially for a large grid of mod- 
els. We calculate the excitation at a number of radial 
points (~ 100 — 120), each of which is effectively treated 
as a separate ID plane-parallel slab. The escape proba- 
bility is therefore given by: 



1 — cxp(— 3t) 
37 ! 



(1) 



where r is the optical depth in the direction of the con- 
sidered ray. The resulting vertical excitation structures 
are subsequently combined and regridded to polar coor- 
dinates, so that it can serve as input for the next step: 
raytracing with RADLite. 

(iii) RADLite (des cribed in a companion paper, Pon- 
toppidan et al. l2009f ) is used to obtain the line spectra. 
It is a raytracer for axisymmetric geometries specifically 
set up to handle the large velocity gradients and small 
turbulent widths which are expected to be present in the 
planet-forming region of circumstellar disks. 

To demonstrate how the model can be used to inves- 
tigate the influence of different water abundance struc- 
tures, two additional physical processes that determine 
the distribution of water vapor are included and de- 
scribed in the following sections. 

4.1. Reduction of the water abundance due to freeze-out 
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Fig. 2. — Dust size distribution with minimum grainsize a m i n = 3 X 10 (im and maximum grainsize a ma x = 40 fim and opacities of 
the silicate (solid) and carbonaceous grains (dashed). Most of the silicate grain mass is in the larger grains (R = 0.5 — 20/^m), while the 
carbonaceous grains have substantial reservoirs in both in small (R < 0.01/xm) and large grains (R > 1/im). 



Freeze-out of water onto grain surfaces can significantly 
reduce the amount of water in the gas phase, and will 
reduce the integrated line emissi on of the lower excita- 
tion lines dMeiierink et all I2008D . The approach from 
iPontoppidanl ( 2006f ) is used to determine the pressure 
and dust temperature where the transition from water 
vapor to ice occurs on grains. The rate of ice mantle 
build-up is given by: 



so high that freeze-out is not important. The freeze-out 
temperature is between T = 100 — 110 K for dense ISM 
cloud conditions, n = 10 4 — 10 6 cm~ 3 , and rises to tem- 
peratures as high as T = 150 — 180 K for midplane den- 
sities in the inner region, R = 0.1 — 0.2 AU, of the disk, 
n = 10 12 — 10 14 cm -3 . If a region is below the freeze-out 
temperature, we adopt a residual gas phase abundance 
of xh 2 o = 10~ 10 to simulate cosmic ray desorption. 



- Rads ~ Rdes, (2) 

where R des = v Q exp(—E/kT dust ) x n H2 cucc x 7 is the 
thermal desorption rate and R a ds — "-h 2 o x iT'd.ust x 
nd 2 x yj3kT gas /mn 2 o x /. The surface binding energy 
of E = 577 3 K is appropriate for water on a crystalline 
ice surface ()Fraser et al.ll200H) . The desorption and ad- 
sorption rates depend on the dust and gas temperature, 
respectively, but freeze-out is only important where the 
gas is coupled to dust, i.e., T gas = Td us t- Non-thermal 
desorption mechanisms are not included. 7 is a factor 
that takes into account the observation (in the labora- 
tory) that H2O only desorbs from the top monolayer (i.e., 
first order desorption for sub-monolayer coverage and ze- 
roth order desorption for multilayers), d — 1 /im is the 
average adopted grains radius, which is representative 
for that used by the radiative transfer models, although 
the precise choice does not affect the results, / = 1 is 
the sticking coefficient, and nd us t and nn 2 are dust and 
water number densities. The pre-exponential factor vq 
is expressed in s" 1 and is roughly the frequency of the 
vibrational stretching mode (A ~3.1 /im). 

In Fig. [31 the density-dependent freeze-out tempera- 
ture is shown. The system is assumed to have reached 
equilibrium, i.e., dni ce /dt = - a reasonable assump- 
tion because the freeze-out timescales are short at the 
high densities in the disk surface (n > 10 8 cm" 3 ). At 
lower densities, n = 10 — 10 6 cm -3 , where equilibrium 
is not reached on time scales shorter than 10 4-5 years and 
where photo-desorption occurs due to direct irradiation 
from the central star, dust temperatures are generally 




100 150 200 250 

Temperature [K] 

Fig. 3. — The density dependence of the freeze-out temperature 
of H2O. Water is mostly frozen out for a give temperature and 
density left to the curve, and mostly in the gas phase when right 
to the curve 



4.2. Reduction of the H2 O abundance in the warm 
atmosphere 

From static chemi s try m o dels (iGlassgold et al.l 12004 
iKamp fc Dullemondl 12004 IWoitke et all l2009f ) it is 
known that the gas and dust are kinetically uncoupled 
in the warm atmosphere of the disk where the ambi- 
ent densities are in the range fin = 10 4 — 10 9 cm~ 3 , 
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Fig. 4. — Left: The e-folding timescales for freeze out are indicated by the black curves for selected scale heights in the disk. The 
Keplerian dynamical timescale is indicated by the dashed line for comparison. Middle: Gas total hydrogen number density. The snow-line 
is indicated as well. Right: Dust temperature. 

leading to local gas temperatures that are significantly 
higher than that of the dust. The gas is at high tem- 
peratures (T ~ 5000 K) in the uppermost unshielded 
regions of the disk and drops to lower temperatures 
(T ~ 100 — 200 K) in deeper layers where the gas is 
coupled to the dust. If this transition zone contains wa- 
ter, it will have a strong effect on the line strengths of 
especially the high excitation lines. The gas/dust decou- 
pling is associated with a reduced water abundance in 
the unshielded part of the disk due to photo-dissociation 
and ion-molecule reactions. A simple parametrized de- 
scription is adopted here using the water abundance 
and temperature st ructure as shown in F ig. 2 at ra- 
dius R ~ 1 AU of iGlassgold et all (|2009f K Warm wa- 
ter is present at temperatures T ~ 300 — 1000 K, but 
reduced by orders of magnitude at very high tempera- 
tures T > 2000 K. The gas temperature and water abun- 
dance depend solely on the ionization parameter C/nn, 
where £ is the total ionization rate and nu the num- 
ber density of hydrogen. The ionization parameter is 
calculated assuming X-rays are the dominated radiation 
source. We adopt a thermal source with Tx — 1 keV 
and Lx = 10 29 erg/s. The total ionization rate for ev- 
ery point on the g rid is calculated using the absor ption 
cross sections from Mo rrison fc McCammonl fl983), and 
taking into account the attenuating column N-r. For ev- 
ery 1 keV of energy absorbed, there are 1 000/37 ~ 27 
ionizations (see e.g. Glassgold et al. 1 19971 for more de- 
tails). Using th i s pres cription, it is possible to scale the 
IGlassgold et al.l (|2009t) results at R = 1 AU to the entire 
disk given a local ionization parameter, and the result- 
ing gas temperature is shown in Fig. [5] The adopted 
thermal and chemical profile can of course be affected by 
vertical mixing, but such considerations are beyond the 
scope of this paper. 

5. RESULTS 

This section is divided into two different parts. The 
first part highlights the differences in H2O line spectra 
from the LTE and non-LTE calculations for the same 
density and temperature distribution, while Section I5~2l 
shows how the line spectra are affected by processes such 
as freeze-out, surface heating, chemistry and transport. 



5.1. LTE versus non-LTE 



Fig. 5. — Gas temperature when heating of the upper atmosphere 
by X-rays is included (Model 5). 

In LTE the level population depends only on the tem- 
perature of the gas. Conversely, in a non-LTE calculation 
all relevant collisional and radiative excitation processes 
must be explicitly included. One consequence of non- 
LTE is that levels may be subthermally excited when 
the ambient densities are lower than the critical densi- 
ties, or superthermally excited when radiative excitation 
dominates the level populations. When a level is sub- 
thermally populated in a particular region of the disk, it 
has a smaller population than in LTE. The ratio between 
the non-LTE and LTE level populations therefore indi- 
cates whether a particular level is thermalized or not. A 
comparison of the non-LTE and LTE face-on integrated 
column densities shows whether the bulk of the gas at 
a certain radius is subthermally excited in a particular 
level. 

Fig. [6] shows the non-LTE / LTE integrated col- 
umn density ratios for all ortho-H^O levels with energies 
E < 7194 K at three different radii with the same disk 
temperature and density distribution and a constant wa- 
ter abundance xh 2 o = 3 x 10~ 4 . It is immediately clear 
from this figure that the LTE approximation is only valid 
at very small radii (R < 0.2 AU). At larger radii the lev- 
els become sub-thermally populated, especially the high 
excitation E > 2000 K lines. The subthermal decrease 
in column density is larger for levels with a higher ex- 
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Fig. 6. — The non-LTE / LTE integrated column density ratios 
for radii R = 0.11, 0.32, and 1.0 AU, for levels in the ground (black) 
and first (red) vibrational states. 

citation energy. Because both density and temperature 
decrease with radius, this results in a smaller emitting 
radius for lines emitting from high excitation levels. The 
collisional excitation rates C u i (T) decrease strongly with 
temperatures from T = 2000 to 200 K in the disk, in 
some cases by more than two orders of magnitude. This 
leads to a steep increase of the critical densities with ra- 
dius. Radiative excitation in part counteracts the lower 
collisional excitation rates in the regions where the gas 
becomes subthermally populated, but for the models in 
this paper, it does not result in a significant increase in 
the integrated line flux. 

Fig. [7] shows the ratios between the non-LTE and LTE 
velocity integrated line fluxes in the Spitzer-IRS band. 
The points are marked with four different colors depend- 
ing on their upper level energy. The ratios are gener- 
ally less than unity, meaning that an assumption of LTE 
will overestimate line fluxes. The figure illustrates that: 
1) the spread in the ratios is very large, 2) the lower 
the upper level energy E up , the closer the ratio is to 
unity, 3) the overall shapes of the LTE and non-LTE 
line spectra (when viewed as an ensemble over the entire 
mid- infrared range) are very different (see also Fig. [5]), 
and (4) the ratio shows only a shallow trend with wave- 
length. The last point is an indication that the upper 
level energies are evenly distributed with line frequency 
thanks to the asymmetric top nature of the water spec- 
trum. This is particularly convenient from an observer's 
perspective as it will usually be possible to observe lines 
with a wide range of properties even if only a limited 
wavelength range is available. 
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Fig. 7. — Ratios of the non-LTE and LTE velocity integrated line 
fluxes for line with upper energy E < 2900 K (top left), 2900 < 
E < 4300 K (top right), 4300 < E < 5800 K (bottom left), and 
E > 5800 K (bottom right) in the Spitzer-IRS band. A distinction 
is made between the ground vibrational level (black crosses) and 
first vibrational level (red triangles). 

Fig. [5] shows the Einstein A coefficients versus the 
non-LTE/LTE line flux ratios, again separated into four 
groups as in Fig. [JJ The figure illustrates the following: 

(i) The line flux ratios decrease for larger Einstein A 
coefficients. This occurs because the critical densi- 
ties are higher for transitions with larger Einstein 
A coefficients (the de-excitation rates are compara- 
ble when the upper level energies are similar). Such 
transitions are therefore farthest from LTE. 

(ii) For this particular model, lines with Einstein A 



coefficients below ~ 10 



are close to LTE. 



(iii) The ratios for transitions within the ground are 
higher than the first vibrational level. This is an in- 
dication that collisions dominate the excitation. If 
fluorescent excitation through the first vibrational 
bending mode were dominant then the line fluxes 
(at the same energy) become significantly larger for 
transitions from the first vibrational level that for 
those from the ground state. 

5.2. Toward a fiducial model 

In this section a model that qualitatively reproduces 
the shape of the observed H2O line spectra is con- 
structed. Fig. [5] shows the line spectra for five differ- 
ent models of non-LTE and LTE spectra convolved to 
two different resolutions, A/AA = 600 (500 km/s, Spitzer 
IRS) and 50000 (6 km/s). All models have the same gas 
density distribution (see Fig. 2J), but have an increasing 
level of complexity of the water distribution and excita- 
tion to accommodate the observables. 

Model 1 is the reference model with a gas tempera- 
ture equal to the dust temperature and a constant water 
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Fig. 8. — Einstein A coefficients versus ratios of the non-LTE 
and LTE velocity integrated line fluxes for line with upper energy 
E < 2900 K (top left), 2900 < E < 4300 K (top right), 4300 < 
E < 5800 K (bottom left), and E > 5800 K (bottom right) in 
the Spitzer IRS band. A distinction is made between the ground 
vibrational level (black diamonds) and first vibrational level (red 
triangles). 

abundance (xh 2 o = 3 x 10~ 4 ) throughout the disk, even 
in regions where water should freeze out. In Model 2 
a density dependent freeze-out of water has been added 
(see Fig. [31 and Section l4~Tj) . The main difference in the 
line spectra between Model 1 and 2 is that some of the 
low excitation lines decrease in intensity. The overall line 
spectrum does not change drastically in the Spitzer range 
with freeze-out added (especially in the non-LTE case). 
This is due to the relatively high excitation temperatures 
of the lines from 10 — 36 fim. as well as the fact that these 
transitions are subthermally excited at the location of 
the surface snow line. However, the freeze-out does affect 
lines of lower excitation, located in the spectr al regime 
between 60 and 200/mi (jMeiierink et al.ll2008l ). relevant 
for Herschel-PACS. Although Model 1 and 2 are able to 
produce bright lines of low excitation energy, similarly 
bright high excitation lines are also required by the ob- 
served spectra, inconsistent with these models. 

One way of increasing the strength of the high excita- 
tion lines is to consider the decoupling of the gas temper- 
ature from the dust in the upper part of the atmosphere 
of the disk. The parametrized treatment of this decou- 
pling, as discussed in Section [4721 is added to Model 2, 
producing Model 3. Both the high and low excitation 
lines increase in intensity (between 20-100 percent), but 
the line to continuum ratios do not approach the obser- 
vations, and the low excitation lines are in fact boosted 
more than the high excitation lines. Furthermore, the 
line-to-continuum ratio is still a factor of 3-5 too low 
with respect to representative Spitzer IRS spectra. For 
instance, the prominent 33 /im line complex now has a 
line to continuum ratio of 0.05:1 at Spitzer-IRS resolu- 
tion, compared to observed values of ~0. 1-0.2:1. 



In order to increase the general line-to-continuum ra- 
tio, Model 4 is introduced, which is the same as Model 
3, except that the gas-to-dust ratio has been increased 
by factor of 10 to simulate strong dust settling. The re- 
sult is that higher gas densities contribute near the disk 
surface leading to increased line intensities, coupled with 
a decrease in the dust continuum. This model shows a 
much better agreement with observations, but many low 
excitation lines are now too bright. Furthermore, the line 
spectrum shows a strong increase in peak line strength 
with wavelength, in contradiction with the observed, al- 
most flat, line integrated intensities. 

The reason that the low excitation lines are bright is 
that they are produced over a wide range of radii in the 
disk. Hence, they can be suppressed if the water va- 
por abundance is drastically lowered beyond a certain 
radius. Static chemical models do predict that the wa- 
ter abundance is lowered by 1-2 orders of magnitude at 
temperatures below ~ 300 K, which might be enough to 
reproduce some of the observed line spectra, but due to 
the high optical depths that remain at low temperature, 
the data indicate that higher depletions may be neces- 
sary. Thus, other mechan isms should be investigated. 
I Stevenson fc Lunine] (fl98§) suggested that water vapor 
(and therefore oxygen) is depleted from the inner disk, 
within the midplane snow line, on time scales shorter 
than the life time of the disk. In this mechanism the wa- 
ter vapor is transported from just within the snow line 
(closer to the central star) to just outside the snow line 
by turbulent diffusion, at which point the water freezes 
out. Essentially, the sharp change in the partial pressure 
of water vapor across the snow line results in a net flux 
of water molecules outwards in the disk. However, this 
radial "cold finger" effect is likely balance d by the in- 
ward migr ation of icy bo dies, as discussed in Sal vk et "all 
(|200cl and|Ciesla| ((20091) . 

The radial cold finger effect assumes that the snow 
line is a vertical line in an isothermal disk. However, 
due to the temperature inversion in the disk atmo- 
sphere dGlassgold et al.|[200l iKamp fc Dullemondl[200l 
iWoitke et al. 2009), the snow line is actually a curve that 
becomes almost parallel to the disk surface over a very 
wide range of the disk, typically the entire planet-forming 
region (0.5-20 AU). This can be seen in Fig. 01 and is 
sketched in Fig. [TU1 Hence, near the surface, the cold fin- 
ger effect will operate in the vertical direction; water va- 
por will diffuse via turbulence to regions below the snow 
line where it will freeze out and take part in the general 
settling of dust to the midplane in a process one might 
call the "vertical cold finger" effect. The important dif- 
ference relative to the radial cold finger effect is that in 
the surface water vapor abundance will not be replen- 
ished by the migration (in the case, the relofting) of solids 
if icy dust grains grow to significant size as they settle and 
as seems to be required by the significant millimeter-wave 
fluxes of many disks, partic ularly those th at are referred 
to as transitional disks (see lBrownl [2007). Further, the 
turbulence will be higher in the surface, thus shortening 
the water depletion time scale. Therefore, we speculate 
that such a mechanism is an effective means of depleting 
surface water above the midplane snow line. If true, this 
would be of considerable importance because it would 
mean that the radial distribution of surface water, as ob- 
served with Spitzer-IRS and future mid- to far-infrared 
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Fig. 9. — Comparison of non-LTE (column a) and LTE (column row) line spectra for constant abundance distribution (model 1), then adding density dependent freeze-out (model 2), 
decoupling gas from dust in warm atmosphere (model 3), increasing the gas to dust ratio from 1280 to 12800 (model 4), and reducing the H2O abundance above mid-plane snowline, 
proposed as vertical cold finger effect (model 5), respectively. 
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facilities, may well trace the midplane snow line. More 
detailed modeling is needed to estimate whether the time 
scales for this proposed process are short enough to oper- 
ate efficiently. It can be noted that a vertical cold finger 
effect truncation above the snow line at ~ 1 AU neatly 
explains why the slab models fit so well when assuming 
the same emitting area for all water lines over a wide 
range of excitation energies. In this way, the water emis- 
sion radius derived from slab models may in fact be a 
measure of the location of the midplane snow line. 

In Model 5 we mimic the vertical cold finger effect by 
setting the water abundance to the freeze-out abundance 
(10~ 10 ) above the midplane snow line of the disk, as il- 
lustrated in Fig. [TO] The resulting spectrum shows that 
the line spectrum has significantly flattened, in general 
agreement with the observed water spectra (see also Fig. 




Fig. 10. — Vertical cold finger effect: Water vapour will diffuse to 
regions below the snowline, where it will freeze out and take part 
in the general settling of dust to the midplane. For the fiducial 
model, the inner rim and the mid-plane snowline are located at 
R=0.07 and 0.7 AU, respectively. 

6. DISCUSSION & CONCLUSIONS 

Driven by the recent detection of hundreds of wa- 
ter emission lines in the mid- infrared f rom disks around 
T Tauri stars by ICarr k Naiital (|2008f) and ISalvk et all 
(2008), a non-LTE 2D model has been constructed 
that qualitatively reproduces the observed spectra. The 
model can be adjusted to obtain detailed fits to indi- 
vidual mid-infrared spectra. Our approach is to derive 
an abundance structure for water based on observational 
constraints, rather than use the prediction of existing 
chemistry or transport models. As a demonstration of 
this technique, the action of a vertical cold finger effect 
is proposed to model the observed truncation of the sur- 
face abundance of water vapor. 
The four most important conclusions are: 
1) A non-thermal excitation treatment of water is es- 
sential in determining the water distribution given an 
observed mid-infrared water spectrum. Assuming LTE 
in the calculation of H2O line spectra is generally in- 
valid, especially for high excitation lines. This is best 
illustrated in Fig. [6] and [9] The critical densities for 



exciting the water lines in the 10 to 36 /im region range 
from n cr it = 10 8 — 10 12 cm -2 , and rapidly increase with 
the upper level energy. In Models 1 and 2, the difference 
between the LTE and non-LTE dominant line intensities 
is a factor of 2 to 3. In these models, the gas temperature 
is assumed to be equal to that of the dust, leading to faint 
high excitation lines. In Model 3, the surface gas temper- 
ature is increased using a param etrization based on de - 
tailed models of the gas heating (jGla ssgold et al. 2009). 
The emission of especially the low excitation lines is now 
boosted by a thin layer of warm H2O in the upper part 
of the disk with lower densities, leading to a much larger 
discrepancy between the non-LTE and LTE case due to 
subthcrmal excitation. To boost the line-to-continuum 
ratio even more, the gas-to-dust ratio in Models 4 and 
5 was increased by a factor of 10 with respect to Model 
1 through 3. The low excitation lines in Model 4 now 
become too bright and therefore a speculative vertical 
cold finger effect was added in Model 5. In this case, 
the difference between the LTE and non-LTE spectrum 
is minimized because the regions of the disk that have 
the largest deviation from LTE contain little/no water 
vapor. In summary, beyond a radius R ~ 0.3 AU, levels 
having transitions in the Spitzer IRS wavelength range 
are sub-thermally excited, leading to biased estimates of 
the temperature and density distribution, emitting ar- 
eas and interpretation of high resolution (~ 3 km/s) line 
profiles when using LTE models. 

2) In order to obtain more line emission from high ex- 
citation water lines (E up > 3000 K), we introduced a 
steep gas temperature increase at high altitudes in the 
disk, which is motivated by both observations and mod- 
els of, e.g., the [Nell] 12.81 fun and [OI] 6300 A lines. 
This increases the line strenghts by 20-100 percent, and 
the low excitation lines are boosted most. This is not 
enough, however, to reproduce the observed Spitzer IRS 
line-to-continuum ratios. 

3) The gas-to-dust ratio must be increased by one to 
two orders of magnitude with respect to the canonical 
ISM ratio of ~ 100 — 200 in order to approach the ob- 
served line strengths and line-to-continuum ratios. An 
increase in gas-to-dust ratio lowers the dust opacity, al- 
lowing lines to be formed at higher densities deeper in 
the disk. An increase in the apparent gas-to-dust ratio 
can be explained by grain growth and dust settling. In- 
creasing the gas-to-dust ratio is a more efficient way to 
increase line-to-continuum ratios than the decoupling of 
gas from dust at high altitudes in the disk because the 
former modification emphasizes high density gas. 

4) Chemic al models (e.g., iGlassgold et all 120091 : 

IWoitke et al.1 [2009) show that it is possible to obtain 
high water abundances in the transition zone from hot 
(T ~ 5000 K) to cold gas (T - 100-200 K) and that the 
range along the vertical axis where the maximum abun- 
dance occurs becomes smaller with radius. Below this 
vertical zone, closer to the midplane, water is less effi- 
ciently formed because the neutral-neutral reactions that 
form it (O + H 2 -► OH + H followed by OH + H 2 -> H 2 
+ H) contain activation barriers. Even so, the predicted 
lower limit to the water abundance (a;n 2 o ~ 10~ 6 ), or a 
suppression by 1-2 orders of magnitude, still produces too 
much emission in the optically thick low excitation lines. 
Rather, a sharp drop in the surface water abundance by 
up to 4 orders of magnitude beyond ~ 1 AU better re- 
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produces the observations. This sharp drop is difficult 
to explain with only static chemical models, although 
it is not completely ruled out that additional chemical 
porcesses could effectively reduce the water abundance. 
We speculate a vertical cold finger effect, a process in 
which water vapor is diffusively transported in the ver- 
tical direction to a point below the snow line where the 
water freezes out and takes part in the general growth 
and settling of solids to the midplane, may provide a vi- 
able explanation (Fig. [T0|) . If this effect does occur and 
is efficient, the radial distribution of water line emission 
from the disk surface can become a direct tracer of the 
location of the midplane snow line. 

A chemical suppression of the water abundance, as 
compared a vertical cold finger process would give rise 
to two completely different disk surface chemistries as 
a function of radius that can be tested observationally: 
If vertical transport is not important and chemical pro- 
cesses along determine the surface water abundance, the 
chemistry is oxygen dominated. Conversely, a verti- 
cal cold finger effect will dramatically reduce the total 
oxygen abundance in the surface, giving rise a carbon- 
dominated chemistry resulting in very different radial 
distributions of organics and other species as observed 
in mid-infrared surface molecular tracers. To test this, 
sensitive high resolution mid-infrared spectroscopy is re- 
quired, such as will be offered by the James Webb Space 
Telescope ( JWST) and the next generation of Extremely 
Large Telescopes (ELTs) , combined with detailed excita- 



tion modeling. 

The present study has qualitatively determined the wa- 
ter temperature and abundance distribution required to 
explain observed mid-infrared spectra of T Tauri disks. 
However, several important aspects that may contribute 
to the formation and appearance of mid-infrared water 
lines from protoplanetary disks are not discussed in this 
paper. 

(i) Intrinsic line broadening: Here we adopted an in- 
trinsic linewidth that is dominated by thermal broaden- 
ing. It is possible that additional broadening from MRI 
driven turbulence would be able to double the intrinsic 
linewidth, although the small predicted values for turbu- 
lent viscosities would indicate that intrinsic line widths 
are indeed dominated by thermal broadening. An in- 
crease of the intrinsic line broadening would reduce the 
opacities in the lines, thus increasing the effective crit- 
ical densities, leading to lowered intensities of the high 
excitation lines. 

(ii) Stellar properties: The central stars in the observed 
Spitzer sample have masses ranging M = 0.3 to 3.0 M©, 
and ages ranging from 1 to 10 Myrs, corresponding to 
luminosities and effective temperatures ranging from L = 
0.1 to 20 L and T cS = 2000- 10000 K, respectively. The 
modeling presented here has not explored dependencies 
on the properties of the central star. 

(iii) Inner holes: The maximum dust temperature will 
decrease when the inner rim is located at larger radii, 
meaning that the high excitation lines will become less 
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prominent in disks with inner holes, while some lower 
excitation lines may actually become brighter (see Pon- 
toppidan et al. 2009). 

(iv) The model setup can easily be extended to other 
molecular species. Work is underway to treat other ob- 
served molecular tracers in the mid-infrared, including 
CO, HCN, C 2 H 2 , OH, etc. 

These aspects will be discussed in a subsequent pa- 
per dedicated to a full parameter study (Meijerink et al. 
2009, in prep.). 
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